{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9978c6fe",
   "metadata": {},
   "source": [
    "# Tutorial 04: computation of the inf-sup constant for the weak imposition of Dirichlet BCs by a Lagrange multiplier\n",
    "\n",
    "In this tutorial we compute the inf-sup constant of the saddle point problem resulting from the discretization of the following Laplace problem\n",
    "$$\\begin{cases}\n",
    "-\\Delta u = f, & \\text{in } \\Omega,\\\\\n",
    " u   = 0, & \\text{on } \\Gamma = \\partial\\Omega,\n",
    "\\end{cases}$$\n",
    "\n",
    "where $\\Omega$ is a ball in 2D, and for which the non-homogeneous Dirichlet boundary conditions are imposed by a Lagrange multiplier.\n",
    "\n",
    "The resulting eigenvalue problem is\n",
    "\\begin{align*}\n",
    "&\\text{find } \\eta, u, \\lambda \\in \\mathbb{R} \\times V \\times M \\text{ s.t. }\\\\\n",
    "&\\begin{cases}\n",
    "\\int_\\Omega \\nabla u \\cdot \\nabla v + \\int_\\Gamma \\lambda v = 0, & \\forall v \\in V,\\\\\n",
    "\\int_\\Gamma u \\mu = \\eta \\int_\\Gamma \\lambda \\mu, & \\forall \\mu \\in M\n",
    "\\end{cases}\n",
    "\\end{align*}\n",
    "where\n",
    "$$\n",
    "V = H^1(\\Omega),\\\\\n",
    "M = L^{2}(\\Gamma).\\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8b80bc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfinx.fem\n",
    "import dolfinx.io\n",
    "import gmsh\n",
    "import mpi4py.MPI\n",
    "import numpy as np\n",
    "import petsc4py.PETSc\n",
    "import slepc4py.SLEPc\n",
    "import ufl\n",
    "import viskex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88b3611f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiphenicsx.fem\n",
    "import multiphenicsx.fem.petsc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65147329",
   "metadata": {},
   "source": [
    "### Geometrical parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9167b8c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 3\n",
    "mesh_size = 1. / 4."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5f95912",
   "metadata": {},
   "source": [
    "### Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0ada126",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.initialize()\n",
    "gmsh.model.add(\"mesh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0970e93c",
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)\n",
    "p1 = gmsh.model.geo.addPoint(0.0, +r, 0.0, mesh_size)\n",
    "p2 = gmsh.model.geo.addPoint(0.0, -r, 0.0, mesh_size)\n",
    "c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)\n",
    "c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)\n",
    "boundary = gmsh.model.geo.addCurveLoop([c0, c1])\n",
    "domain = gmsh.model.geo.addPlaneSurface([boundary])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef2f0599",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.model.geo.synchronize()\n",
    "gmsh.model.addPhysicalGroup(1, [c0, c1], 1)\n",
    "gmsh.model.addPhysicalGroup(2, [boundary], 0)\n",
    "gmsh.model.mesh.generate(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f26829f",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(\n",
    "    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bf09b8a-b87b-4fc3-966f-beed80d3a426",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create connectivities required by the rest of the code\n",
    "mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d06090a",
   "metadata": {},
   "outputs": [],
   "source": [
    "facets_Gamma = boundaries.indices[boundaries.values == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d5e7bb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07a9fc20",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh_tags(mesh, boundaries, \"boundaries\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e721275b",
   "metadata": {},
   "source": [
    "### Eigenvalue problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f3e7f17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function space\n",
    "V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2cc8c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define restrictions.\n",
    "dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)\n",
    "dofs_V_Gamma = dolfinx.fem.locate_dofs_topological(V, boundaries.dim, facets_Gamma)\n",
    "restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, dofs_V)\n",
    "restriction_V_Gamma = multiphenicsx.fem.DofMapRestriction(V.dofmap, dofs_V_Gamma)\n",
    "restriction = [restriction_V, restriction_V_Gamma]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0073a5e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define trial and test functions\n",
    "(u, l) = (ufl.TrialFunction(V), ufl.TrialFunction(V))\n",
    "(v, m) = (ufl.TestFunction(V), ufl.TestFunction(V))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "236b1f94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define problem block forms\n",
    "a = [[ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, - ufl.inner(l, v) * ufl.ds],\n",
    "     [- ufl.inner(u, m) * ufl.ds, None]]\n",
    "b = [[None, None],\n",
    "     [None, - ufl.inner(l, m) * ufl.ds]]\n",
    "b[0][0] = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(u, v) * ufl.dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "291922b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assemble lhs and rhs matrices\n",
    "A = multiphenicsx.fem.petsc.assemble_matrix_block(\n",
    "    dolfinx.fem.form(a), bcs=[], restriction=(restriction, restriction))\n",
    "A.assemble()\n",
    "B = multiphenicsx.fem.petsc.assemble_matrix_block(\n",
    "    dolfinx.fem.form(b), bcs=[], restriction=(restriction, restriction))\n",
    "B.assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21f9dc48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve\n",
    "eps = slepc4py.SLEPc.EPS().create(mesh.comm)\n",
    "eps.setOperators(A, B)\n",
    "eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)\n",
    "eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)\n",
    "eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.TARGET_REAL)\n",
    "eps.setTarget(1.e-5)\n",
    "eps.getST().setType(slepc4py.SLEPc.ST.Type.SINVERT)\n",
    "eps.getST().getKSP().setType(\"preonly\")\n",
    "eps.getST().getKSP().getPC().setType(\"lu\")\n",
    "eps.getST().getKSP().getPC().setFactorSolverType(\"mumps\")\n",
    "eps.solve()\n",
    "assert eps.getConverged() >= 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60604a4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract leading eigenvalue\n",
    "eigv = eps.getEigenvalue(0)\n",
    "r, i = eigv.real, eigv.imag\n",
    "assert abs(i) < 1.e-10\n",
    "assert r > 0., \"r = \" + str(r) + \" is not positive\"\n",
    "print(\"Inf-sup constant: \", np.sqrt(r))\n",
    "assert np.isclose(np.sqrt(r), 0.125429, rtol=np.inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b82f3acd",
   "metadata": {},
   "outputs": [],
   "source": [
    "eps.destroy()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
